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We present a calculation of the differential two jet cross section in e + e~ annihilation through 
next-to-next-to- leading order in the strong coupling constant a s - The calculation is performed 
using a new method for dealing with real radiation suggested recently by us in Q. For the first 
time, the two jet event rate is computed directly, without any reference to the inclusive cross-section 
e + e~ — > hadrons. We also calculate the energy distribution of the leading jet in e + e~ — » 2 jets and 
find significant modifications of the shape of this distribution at NNLO. 



PACS numbers: 

High-energy physics will begin to explore a new energy 
frontier when the Large Hadron Collider at CERN turns 
on in 2007. Our understanding of physics at very small 
distances will dramatically improve. However, a detailed 
investigation of the new physics we discover will require a 
careful study of Standard Model backgrounds, detector 
responses, and other similar issues. Since many layers 
separate interesting physics from raw experimental data, 
a dedicated effort is required to fully utilize LHC results. 
There have been significant advances towards this goal in 
the past few years; we now have an increased understand- 
ing of parton distribution functions and jet algorithms, 
improved Monte Carlo event generators, methods for au- 
tomating next-to-leading order calculations with large 
number of external legs and, finally, new technology for 
next-to-next-to-leading order (NNLO) computations. 

NNLO calculations are certainly not required for all 
processes at the LHC or existing colliders; however, there 
are a few situations in which NNLO calculations are 
highly desirable. These include processes for which the 
one-loop corrections are abnormally large (e.g. the pro- 
duction of the SM Higgs boson at hadron colliders 0) 
or for measurements in which high experimental preci- 
sion is either achieved (e.g. the a s determination from 
the three-jet event rate in e + e~ annihilation Q or the 
W mass measurement at the Tevatron) or expected (e.g. 
W and Z boson production at the LHC 3). These cal- 
culations should also inform us how accurate NLO cal- 
culations really are, beyond the standard checks of sta- 
bility with respect to renormalization and factorization 



scales variations. We should learn to estimate the signifi- 
cance of NNLO corrections without performing the calcu- 
lations, given the required precision for an observable and 
the kinematic regions in which it is measured. For this 
purpose, exclusive NNLO calculations are needed, since 
experimental cuts on the final state can have a strong im- 
pact on the convergence of perturbative expansion. Un- 
fortunately, not a single calculation of a fully differential 
QCD observable at NNLO has been performed, either for 
lepton or hadron colliders. 

In this Letter we remedy this situation and present 
the calculation of the two jet cross-section in e + e~ an- 
nihilation at NNLO in perturbative QCD. Although jets 
and their properties have been studied very extensively 
at lepton colliders pg , we believe that such calculation is 
important for the following reasons: 1) it is the first-ever 
calculation of a fully differential observable at NNLO; 
2) although the total rate for two jet events in e + e~ 
annihilation is known through NNLO from indirect cal- 
culations, our results for distributions in two jet events 
are new; 3) this calculation is possible because of a new 
method we recently suggested for handling real radiation 
in hard processes 0]; it is important to demonstrate its 
efficiency by applying it to a non-trivial example. 

There is a strong correlation between the complexity 
of higher order calculations and the level of exclusiveness 
desired. Traditionally, it was thought that calculations 
at higher orders are difficult because of multi-loop inte- 
grals. A significant effort therefore went into developing 
flexible, easily automated methods for performing higher 
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loop computations j6|. As a result, calculations up to 
two loops are no longer prohibitively difficult; for exam- 
ple, the two-loop virtual corrections for 1 — > 3 and all 
partonic 2 — > 2 processes at hadron colliders have been 
computed |7j. Surprisingly, the major obstacle in obtain- 
ing differential results at higher orders are tree-level pro- 
cesses with additional final-state partons. While it is very 
easy to write down the corresponding matrix element, it 
is not possible to integrate it numerically over the re- 
stricted (exclusive) phase-space without first extracting 
the singular structure of the integrand in the soft and 
collinear limits. Analytic integrations also become ex- 
tremely difficult because of the arbitrariness of final-state 
cuts. This problem has been successfully solved at NLO 
using both the slicing and dipole subtraction methods 
8, 9]. Attempts have been made recently to generalize 
the dipole formalism to NNLO ^(j ; so f ar i they have not 
completely succeeded. It is therefore productive to look 
for an alternative method of dealing with multi-particle 
final states in the presence of arbitrary constraints on 
their phase space. 

What are the ideal features of such a method? Given 
the complexity of higher order calculations, it should 
satisfy the following requirements: 1) the singularities 
should be extracted in an algorithmic fashion; 2) the 
method should be easy to automate; 3) it should be easily 
generalizablc, at least in principle, to arbitrary numbers 
of partons in the final state; 4) the method should work 
efficiently in the presence of arbitrary constraints on the 
final state; 5) it should lead to a fast and accurate nu- 
merical evaluation of physical quantities. 

We have proposed such a method recently in pj. We 
now briefly describe its salient features. Consider a per- 
turbative tree-amplitude M with n particles in the final 
state. Its contribution to the differential cross-section 
can be written as 

do-W = J dT n \M\ 2 J({ Pi }), (1) 

where dF„ denotes the n-particle phase-space and 
J({pi}) imposes restrictions on the final state (e.g. the 
jet algorithm) which define the experimentally observed 
process. Throughout this Letter we use dimensional reg- 
ularization (with d = 4— 2e dimensions) for both infrared 
and collinear divergences. If all particles are well sepa- 
rated (resolved), \M\ 2 is finite; however, when the inte- 
gration in Eq.JTJ is attempted, there are divergences as- 
sociated with soft and collinear kinematic configurations. 
A direct numerical integration of Eq. Q is therefore not 
possible. 

In Ref. Q we demonstrated that by mapping the in- 
variant masses onto algorithmically chosen new integra- 
tion variables, it is possible to extract the e poles explic- 
itly. We then obtain an expansion in e, 

2(n-2) 

d- (n) = E p^FS + 0(e), ( 2 ) 

where the coefficients dFk are well-defined, e-independent 
multi-dimensional integrals for a generic function J . It 



needs to be specified only at the stage of numerical eval- 
uation. This allows us to derive results for arbitrary jet 
algorithms and experimental observables. 

It is relatively easy to derive such an expansion at 
NLO, where at most one parton can become unresolved. 
In this case, trivial mappings of the phase-space vari- 
ables onto variables with range from to 1 produce in- 
tegrals with the following singular structure: 

l 

h= J dxdy z e -y- x J{x,y). (3) 
o 

All singularities in the above integral can be extracted 
by writing x c_1 = 8{x)/e + [l/x} + + e[ln(x)/x]+ + for 
both x and y in the integrand. 

Beyond NLO, two or more partons may become unre- 
solved, which gives rise to a more complicated structure 
of overlapping singularities. Typically, we find integrals 
similar to 

l 

h= J dxdy j^-^ J(x,y). (4) 
o 

The procedure described above does not work because 
the singularities are not factorized. To solve this prob- 
lem, we apply the technique of sector decomposition ; 
we divide the integration region in Eq.Q) into patches 
with a definite ordering of the integration variables (x < 
y and y < x) and reweight all variables in each patch so 
that the integrations again range from to 1 . This leads 
to factorization of the singular limits. This procedure 
can be completely automated. The same method should, 
in principle, work for any number of particles in the final 
state, both massless and massive, and for any restrictions 
on the final-state phase space. 

For e + e~ — > 2 jets through NNLO, the largest multi- 
plicity of particles in the final state is four (e.g. e + e~ — > 
qqgg). A parameterization of the 1 — > 4 particle phase- 
space in terms of five independent variables which is suit- 
able for sector decomposition and extracting infrared di- 
vergences was given in [J. In the same reference, we 
gave a more detailed description of the method, and con- 
sidered a number of relatively simple examples. In this 
Letter we apply the method to a fully realistic and non- 
trivial problem - the calculation of the e + e~ — > 2 jets 
cross section at NNLO. Traditionally, the inclusive 2-jet 
rate is calculated at NNLO indirectly, by first comput- 
ing the total inclusive cross-section for e + e~ — > hadrons 
and then subtracting from it the e + e~ — > 4 jets and 
e + e~ — > 3 jets cross-sections at LO and NLO, respec- 
tively. This paper presents the first direct calculation of 
the 2-jet rate at NNLO. Using our method, we can also 
obtain differential results at NNLO, which can not be 
derived indirectly. We illustrate this by computing the 
energy distribution of the leading jet in e + e~ — ► 2 jets 
through NNLO. 

The cross-section for e + e~ annihilation into hadrons 
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through order 0(a 2 s ) can be written as 



where <jq = 47raQED J2 q Qq/ S i s the tree level cross- 
section for e + e~ — » qq, y/s is the center of mass energy, 
a s = a s (s) is the MS QCD coupling constant and the co- 

efficients Cf describe jet cross-sections in various orders 
in perturbation theory, as indicated by the Kronecker 
symbols. From inclusive calculations of the cross-section 
we find 
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1.99- 0.1157V/, 



where Nf is the number of massless fermion flavors. 
An important goal of this Letter is the calculation of 

(2) 

the coefficient C 2 , the NNLO correction to the two- 
jet production rate. For this, we need the two-loop vir- 
tual correction to e + e~ — > qq, the one- loop correction 
to the e + e _ — > qqg process, and the tree level processes 
e + e~ —>■ qqgg and e + e~ — * qqq\q\. We also require the 
coupling constant renormalization of the NLO result. At 
order 0{a 2 s ), all of these processes contain divergent con- 
tributions to the two jet cross-sections; the highest sin- 
gularity is 1/e 4 . The singularities cancel when individual 
contributions are combined to form a physical observable. 

The two-loop virtual corrections to e + e~ — > qq are 
well-known [13. We have outlined above how the tree- 
level four parton final state is handled in our approach. 
We note that a global parameterization of the four par- 
ticle phase space, which we used in pj for the TV/ terms, 
leads to large analytic expressions which are difficult to 
evaluate numerically. We found it much more conve- 
nient to select a different parameterization for the in- 
variant masses in each individual term, thereby reducing 
the number of sector decompositions required to extract 
the singularities. This choice of parameterization can be 
done automatically once the basis topologies appearing 
in the matrix element are identified. With this clever 
choice of parameterization, the size of the computer code 



for the fully differential NNLO e H 



2 jets process 



is not much larger than what we have found in simpler 
examples in pj. The required CPU time is also not very 
large; to achieve the precision on the jet rates presented 
in this paper, about four hours are needed on a PC with 
a 3 GHz Pentium 4 processor. 

We now briefly comment on the calculation of the one- 
loop corrections to the qqg final state. It might seem that 
a different technique is needed to handle this contribu- 
tion, since a virtual loop integration is involved. How- 
ever, this is not the case 0. Once the virtual loop in- 
tegrals are expressed through Feynman parameters, they 
can be treated identically to phase-space integrals. We 



found it convenient to express them through standard 
hypergcometric functions, and use the one-dimensional 
integral representation for the hypergeometric functions, 
together with the three-parton phase space parameteri- 
zation. Although this procedure is not necessary, it is 
useful because it provides an economical input for sector 
decomposition. 

Since our approach to the problem is numerical, in- 
cluding the cancellation of 1/e poles, we must consider 
issues of numerical accuracy. The simplest check is the 
comparison of the direct and indirect evaluations of the 
total two-jet event rate. The indirect result is obtained 
by taking the difference between the 0(a s s ) contribution 
to the inclusive cross-section, given in Eq. JHJ, and sub- 
tracting from it the four-jet cross-section at LO and the 
three-jet cross-section at NLO. Both of these quantities 
are computed in our code. We use the JADE algorithm 
to identify jets in the final state. However, the jet 
definition is an independent subroutine in our code that 
can be trivially changed if desired. Choosing the jet sep- 
aration parameter for the JADE algorithm y cut = 0.1, 
we obtain 
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(-49.2 ± 0.4) + (1.7974 ± 0.0011) A/, (7) 



where the errors denote our integration uncertainties for 
the 3 and 4 jet cross sections. A direct computation of 
the same quantity yields 
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We have included the integration errors found during an 
actual run for the 1/e poles to demonstrate the level 
of cancellation; the magnitudes indicated for the higher 
poles are typical of results found using our code. Com- 
paring Eq.Q and Eq.JSJ), we conclude that our program 
provides a precision on the finite part of the NNLO cor- 
rection to the two-jet rate better than 1%. We also con- 
clude that our numerical cancellation of 1/e poles works 
very efficiently. These features do not change signifi- 
cantly when the jet separation parameter y cut is varied. 

Our approach permits us to also compute differential 
distributions in addition to the total rate. As an ex- 
ample, we present below the energy distribution of the 
leading jet in two jet events at NNLO. At leading order, 
this distribution is simple; since two massless quarks are 
produced, each jet contains half of the total energy. The 
distribution becomes more interesting at NLO, when it 
becomes possible for one of the jets to have an invari- 
ant mass different from zero. At NNLO, configurations 
when the invariant masses of both jets are different from 
zero appear for the first time. We compute this distribu- 
tion by a simple modification of the jet function; after an 
event is identified as a two jet event, the energies of the 
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Ycut = 0.1 



LO 



FIG. 1: Bin-integrated energy distribution for y cut =0.1. The 
fractions of events in each energy bin are shown. The dotted 
histogram denotes the LO result, the dashed histogram the 
NLO result, and the solid histogram the NNLO result. 



two jets are computed and the jet with the largest energy 
is identified. This number is then stored in the appropri- 
ate bin of a histogram. The corresponding bin-integrated 
energy distribution is shown in Fig. 1 for y cut = 0.1, 
Nf = 5, and a s = 0.121. The distribution is significantly 
distorted by NNLO QCD corrections; the corrections are 
large for this y cu t, and many situations that look like 
three and four jet events are identified as two jet events. 
Smaller y cut choices lead to large logarithms in the per- 
turbative expansion that invalidate the fixed order result. 



In conclusion, we have presented the first calculation 
of the NNLO corrections to a fully differential observ- 
able in QCD. We have demonstrated our approach using 
the non-trivial example of e + e~ — > 2 jets. We have com- 
puted the NNLO corrections to the energy distribution of 
the two jets in e + e~ annihilation, and have shown that 
the shape of the distribution changes when the NNLO 
corrections are included. Our method allows other phe- 
nomenologically interesting distributions in 2-jet events 
to be easily computed; these will be discussed elsewhere. 
Since our approach is fully numerical, we have presented 
convincing evidence that reasonable precision and control 
of numerical stability can be achieved. The method we 
developed for this calculation is quite flexible; it general- 
izes straightforwardly to an arbitrary number of partons 
in the final state, both massive and massless. Given un- 
limited computing resources, it provides a complete so- 
lution to the problem of real radiation at higher orders 
in perturbative QCD. In practice, significant effort and 
some ingenuity will be required to apply it to more com- 
plicated processes of direct phenomenological relevance. 
We look forward to this challenge. 
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